NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CALIFORNIA 


THESIS 


Approved  for  public  release;  distribution  unlimited. 

19910121  205 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 


PubUc  reporting  burden  for  this  coUection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction,  searching  existing  data 
sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  coUection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other 
aspect  of  this  coUection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and 
Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Oftice  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188) 
Washington  DC  20503.  _ _ _ _ — 


AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE 

September  1996 

TITLE  AND  SUBTITLE  Single  Source  Error  Ellipse  Combination 
AUTHOR(S)  Joseph  R.  Orechovesky 


PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 
Naval  Postgraduate  School 
Monterey  CA  93943-5000 


REPORT  TYPE  AND  DATES  COVERED 
Master's  Thesis 

5.  FUNDING  NUMBERS 


12b.  DISTRIBUTION  CODE 


8.  PERFORMING 
ORGANIZATION 
REPORT  NUMBER 


9  SPONSORING/ MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES)  10.  SPONSOiaNG/ MONITORING 

AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the 

official  policy  or  position  of  the  Department  of  Defense  or  the  U.S.  Government. _ 

12a.  DISTRIBUTION/ AVAILABILITY  STATEMENT  '  12b.  DISTRIBUTION  CODE 

Approved  for  public  release;  distribution  is  xmlimited. 

13.  ABSTRACT  (maximum  200  words) 

There  are  a  number  of  military  applications  in  which  the  geographic  location  of  a  signal  of  interest  is  of 
prime  importance  to  the  ability  of  a  unit  to  fulfill  its  mission.  The  accuracy  of  the  geographic  fix  provided 
to  the  warfighter  can  directly  affect  the  success  or  failure  of  a  mission.  One  method  to  improve  the 
accuracy  of  existing  systems  is  to  use  the  weighted  average  of  a  number  of  intercepts.  Each  intercept  is 
manifested  as  an  error  ellipse  comprised  of  a  latitude,  longitude,  semi-axes,  heading  and  a  related  Chi- 
'  squared  distributed  probability.  Individual  error  ellipses  can  be  Viewed  as  a  quadratic  surface 
perpendicular  to  the  x,y  plane  of  a  bivariate  normal  distribution,  the  z-axis  intersection  of  which 
corresponds  to  a  Chi-squared  value.  By  transforming  the  individual  error  ellipses  to  their  related  location 
^  covariance  matrices,  Gaussian  statistics  may  be  used  to  produce  a  single  location  ellipse  that  combines 
information  from  several  two-dimensional  target  location  ellipses.  By  providing  a  means  to  fuse  data  from 
•  a  given  source  the  warfighter  or  analyst  will  be  able  to  more  accurately  assess  a  threat  and  respond. 


14.  SUBJECT  TERMS  error,  ellipse,  geolocation 


17.  SECURITY  CLASSIFICATION  OF 
REPORT 
Unclassified 


18.  SECURITY  CLASSIFICATION  19,  SECURITY  CLASSIFICA- 
OF  THIS  PAGE  TION  OF  ABSTRACT 

Unclassified  Unclassified 


15.  NUMBER  OF 
PAGES  82 

16.  PRICE  CODE 

"20  UMITAHON  OF 
ABSTRACT 
UL 


1 


11 


Approved  for  public  release;  distribution  is  unlimited. 


SINGLE  SOURCE  ERROR 
ELLPSE  COMBINATION 

Joseph  R.  Orechovesky  Jr. 

Lieutenant,  United  States  Navy 
B.S.,  University  of  South  Carolina,  1989 

Submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  SYSTEMS  TECHNOLOGY 
(Space  Systems  Operations) 

from  the 


Author; 


Approved  by  ; 


NAVAL  POSTGRADUATE  SCHOOL 
September  1996 


/  TrkOAnli  T?  Orprlinvpckv  Tr 


Joseph  R  Orechovesky  Jr. 


Vicente  Garcia,  Thesis  Advisor 


Gerry  Baumgartner,  Secured  Reader 


Rudolf  Panholzer,  Chairman,  SpaVcwstems  Academic  Group 


iii 


IV 


ABSTRACT 


There  are  a  number  of  militar7  applications  in  which  the  geographic  location  of  a 
signal  of  interest  is  of  prime  importance  to  the  ability  of  a  unit  to  fulfill  its  mission.  The 
accuracy  of  the  geographic  fix  provided  to  the  warfi^ter  can  directly  affect  the  success  or 
failure  of  a  mission.  One  method  to  improve  the  accuracy  of  existing  systems  is  to  use  the 
weighted  average  of  a  number  of  intercepts.  Each  intercept  is  manifested  as  an  error  ellipse 
comprised  of  a  latitude,  longitude,  semi-axes,  heading  and  a  related  Chi-squared  distributed 
probability.  Individual  error  ellipses  can  be  viewed  as  a  quadratic  surface  perpendicular  to 
the  x,y  plane  of  a  bivariate  normal  distribution,  the  z-axis  intersection  of  which 
corresponds  to  a  Chi-squared  value.  By  transforming  the  individual  error  ellipses  to  their 
related  location  covariance  matrices,  Gaussian  statistics  may  be  used  to  produce  a  single 
location  ellipse  that  combines  information  from  several  two-dimensional  target  location 
ellipses.  By  providing  a  means  to  fuse  data  from  a  given  source  the  warfi^ter  or  analyst  will 
be  able  to  more  accurately  assess  a  threat  and  respond. 


VI 


TABLE  OF  CONTENTS 


I.  INTRODUCTION . 1 

A.  JUSTIFICATION . 1 

B.  SCOPE . 1 

C.  CHAPTER  OUTLINE . 2 

II.  STATISTICAL  THEORY . 3 

A.  RELATIONSHIP  OF  THE  ERROR  ELLIPSE  TO  A  BIVARIATE  NORMAL 

DISTRIBUTION . 3 

B.  BAYES’  THEORM . 8 

C.  METHODOLOGY . H 

III.  CONCEPT  DEVELOPMENT . 15 

A.  COVARIANCE  MATRICES . 15 

1.  Scaled  Covariance  Matrix . 15 

2.  2x2  to  3x3  Covariance  Matrix  Conversion . 18 

3.  Selection  of  Eigenvalues  and  Eigenvectors . 20 

a.  Eigenvalues . 20 

b.  Identities . 21 

c.  Eigenvectors . 22 

d.  Eigenvector  Matrix . ! . 31 

4.  Variable  Relationships . 32 

B.  N+1  COMBINATORIAL  EFFECTS . 36 

C.  ELLIPSE  ROTATION . 39 

D.  ELLIPSE  HEADING  RELATIVE  TO  NORTH  POLE . 44 

IV.  COMBINATION  ALGORITHM . 51 

A.  PROCESS  OUTLINE . 51 

B.  PROCESS  FLOW  CHART . 61 

V.  SUMMARY  &  CONCLUSIONS . 63 

A.  AREAS  FOR  FUTURE  STUDY . 63 

B.  CONCLUSION . 65 

LIST  OF  REFERENCES . 67 

INITIAL  DISTRIBUTION  LIST . 69 


vii 


viii 


LIST  OF  FIGURES 


Number  Pcige 

Figure  2.1  COVARIANCE  ELLIPSES . 6 

Figure  2.2  OPTIMAL  ESTIMATE  OF  TARGET  LOCATION . 14 

Figure  3.1  PROJECTION  GEOMETRY . 16 

Figure  3.2  GEOMETRY  OF  SCALING  MATRIX . 17 

Figure  3.3  ELLIPSE  ORIENTATION . 33 

Figure  3.4  COMBINATION  EXAMPLE  (a) . 37 

Figure  3.5  COMBINATION  EXAMPLE  (b) . 37 

Figure  3.6  ROTATION  COORDINATE  SYSTEM . 40 

Figure  3.7  HEADING  ERROR . 45 

Figure  3.8  HEADING  BETWEEN  POINTS  BAND  A . 46 

Figure  3.9  HEADING  BETWEEN  POINTS  A  AND  B . 49 

Figure  4.1  COORDINATE  TRANSFORMATION . 54 

Figure  4.2  NORTH  POLE  ROTATION . 55 

Figure  4.3  PROCESS  FLOW  CHART . 61 


IX 


X 


ACKNOWLEDGMENTS 


I  would  like  to  express  my  sincerest  thanks  to  all  the  people  that  made  the  completion  of 
this  thesis  possible.  First,  I  want  to  express  my  deepest  gratitude  to  Dr.  Gerry  Baumgartner 
for  giving  me  his  confidence,  enduring  patience,  guidance  and  encouragement.  Without  you 
this  work  could  not  have  been  possible  -  thank  you.  A  special  thanks  to  Professor  Vicente 
Garcia,  no  student  could  find  a  more  caring  or  motivating  advisor.  Finally,  I  want  to  say 
thank  you  to  the  many  Naval  Postgraduate  School  Faculty  and  Staff  that  put  on  a  smile  and 
went  out  of  their  way  to  aid  me  in  one  endeavor  or  another. 


XI 


xii 


I.  INTRODUCTION 


A.  JUSTIFICATION 

There  are  a  number  of  military  applications  in  which  the  geographic  location  of  a 
signal  of  interest  is  of  prime  importance  to  a  unit’s  ability  to  fulfill  its  mission.  The  accuracy 
of  the  geographic  fix  provided  to  the  warfi^ter  can  directly  affect  the  success  or  failure  of  a 
mission.  It  is  for  this  reason  we  must  continually  strive  to  improve  our  ability  to  generate 
the  most  accurate  and  condensed  geographic  solution  from  all  available  data. 

This  thesis  describes  an  algorithm  to  produce  a  single  location  ellipse  that  is  the 
combination  of  several  two-dimensional  target  location  ellipses.  The  combination  algorithm 
presented  in  this  document  is  limited  to  a  stationary  single  source  for  data  set  production. 
The  purpose  of  this  study  is  to  assist  the  Naval,  Command,  Control  and  Ocean  Surveillance 
Center  Research,  Development,  Test  and  Evaluation  Division’s  (NRaD),  Integrated  Satellite 
&  Link  Communications  Division  (Code  841)  in  the  documentation  of  ihe  combination 
algorithm.  This  effort  is  in  support  of  the  Classic  Crystal  project  developed  by  NRaD. 

B.  SCOPE 

The  scope  of  this  thesis  involves  the  development  of  the  combination  algorithm.  A 
mathematically  rigorous  presentation  of  assumptions  made  in  the  algorithm  and  of  the 
supporting  theory  will  be  presented.  If  an  appropriate  method  is  available,  verification  of 
developed  mathematical  equations  is  supplied.  It  is  not  within  the  scope  of  this  project  to 
consider  the  effect  of  system  dependent  biases  (e.g.  atmospheric  delays,  refiraction)  on  error 
ellipses  produced  by  different  systems.  For  this  reason  the  data  for  this  project  is  limited  to 
error  ellipses  produced  by  a  single  source. 
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C.  CHAPTER  OUTLINE 

Chapter  II  presents  the  underlying  statistical  theory  upon  which  the  Combination 
Algorithm  is  based.  The  chapter  culminates  with  the  formula  for  the  optimal  estimate  of 
the  target  location  as  presented  by  J.A.  Roecker,  who  based  his  work  on  that  of  Nelson 
Blachman.  The  concepts  of  a  bivariate  normal  distribution  and  application  of  Bayes’ 
Theorem  are  covered  as  they  directly  support  the  formulation  of  the  optimal  estimate. 
Chapter  III  concerns  itself  with  a  number  of  topics  all  related  to  applying  the  optimal 
estimate  of  the  target  location  developed  in  Chapter  II  to  a  real  world  scenario. 
Development  of  the  covariance  matrices,  the  effects  of  adding  an  observation  to  an  existing 
solution,  ellipse  rotation  and  heading  issues  are  covered  in  rigorous  detail.  In  Chapter  FV 
the  theory  and  concepts  of  the  previous  two  chapters  are  put  to  use  in  the  presentation  of 
the  Combination  Algorithm.  The  process  is  detailed  in  a  step  by  step  manner  that  can 
easily  be  applied  to  a  number  of  platforms.  Finally,  Chapter  V  offers  conclusions  and 
recommendations  for  further  study. 
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II.  STATISTICAL  THEORY 


A.  RELATIONSHIP  OF  THE  ERROR  ELLIPSE  TO  A  BIVARIATE 

NORMAL  DISTRIBUTION 


Each  of  the  target  location  error  ellipses  to  be  combined  to  produce  a  single  error 
ellipse  will  be  defined  by  the  following  parameters;  center  position,  heading  and  magnitude 
of  semi-axes  and  probability  level.  These  ellipses  represent  an  area  that  has  the  specified 
probability  that  the  target  lie  within  it’s  bounds.  We  will  begb  by  defining  a  joint  normal 


distribution  of  two  variables  x  = 


as 


^(x)  =  k*  expj-  ^(x  -  p)’^b(x  -  n) 


(2.1) 


where  k  is  determined  by  the  normalization  of  the  total  probability  to  one,  and  x  are  2- 


component  Column  vectors  and  B  a  (2x2)-matrix.  The  vector  \i  — 
of  the  center  of  symmetry  of  the  distribution 


yMi) 


gives  the  location 


J|(x-p)<!Kx>&A2  =0 


(2.2) 


The  expected  or  mean  value  of  a  random  vector  x  is 


E(l)=  J  = 


(2.3) 


The  covariance  matrix  of  the  distribution  is  defined  as 

C(i)  =  E(x  -  - 1*)’  =  [e(^,  -  Hh  -  ^/)1 
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(Anderson, 1984).  The  i  th  diagonal  element  of  this  matrix,  e(x,  -//,  )  ,  is  the  variance  of 
and  the  ij  th  off-diagonal  element,  e(x,.  -  ^)(xj  -  ,  is  the  covariance  of  x,  and  Xj, 

i^j.  Note  that  since  e(x, -/4){x^. -//,)  =  e(x^. -yt/^)(x; -a)  the  covariance  matrix  is 
symmetric.  With  this  said  the  covariance  matrix  can  be  written  as 


C  =  B  ‘ 


(2.5) 


/  2 

f"  2 

<7,2 

0-1 

^12 

II 

2 

— 

2 

(J2J 

V^12 

a^J 

(2.6) 


where  we  have  used  fhe  symmetry  of  the  covariance  matrix  (i.e.  CTjj  -  cr^j).  By  matrix 
inversion  we  obtain 


B 


1 

0-2 

<^1^2 

V  O’jj 

erf  J 

(2.7) 


To  facilitate  the  next  step  in  our  development  of  the  ellipse  of  covariance  we  will  introduce 
two  reduced  variables 


M.  = 


i  =  1,2 


with  the  property  var(«i)  =  var(«2)  =  1,  and  die  correlation  coefficient 

P  =  -^  =  cov(m,,M2). 

Equation  2.1  now  takes  the  reduced  form 

,2/2)  =  k*  exp(^-  ^«'^B«j 


(2.8) 


(2.9) 


(2.10) 
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with 


B  = 


\-p^ 


1  -p 

u 


(2.11) 


Lines  of  constant  probability  density  are  determined  by  requiring  the  exponent  in  Equation 
2.10  to  be  constant,  as  shown  in  the  following  equation: 


+«:  -  2u,u,p)  =  c 


(2.12) 


If  we  choose  the  constant,  c=l,  then  in  terms  of  the  origbal  variables.  Equation  2.12 
becomes 


n  —  A  i 


(2.13) 


(Brandt,  1983).  This  can  also  be  written  as  a  sum  of  squares  of  two  stochastically 
independent  variables 


1 

(^1  -  a)'  .  -  a  ^2  -  a  ,  (^2  -  a)' 

( 

(  \ 

A 

1-p' 

-  -  ip  +2 

<T,  Oi  0-2  0-2 

+ 

IPvJ 

(2.14) 


which  are  distributed  as  with  two  degrees  of  freedom.  This  is  the  equation  of  an  ellipse 
with  the  center  located  at  ju^,  pz-  The  semi-axes  of  the  ellipse  have  an  an^e  B  with  respect 
to  the  Xi,  Xj  axes.  This  angle,  and  the  magnitude  of  the  semi-major  axis  a  and  the  semi¬ 
minor  axis  bean  be  derived  from  Equation  2.13  using  the  properties  of  conic  sections: 


tan2^  = 


2pc^^cr^ 


(2.15) 


_ _ 

cxl  cos^  B-lpa^  sin0cos^-t-  of  sin^  B 


(2.16) 
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The  ellipse  we  have  just  described  is  known  as  the  ellipse  of  covariance  of  the 
bivariate  normal  distribution  or,  in  the  context  of  this  thesis,  the  error  ellipse.  The  ellipse 
will  always  fall  within  the  rectangle  defined  by  points  (///,  JU;^  and  the  standard  deviations 
(Oi_  It  can  been  seen  from  Figure  2.1  that  the  ellipse  will  touch  the  rectangle  at  four 
points  and  in  the  extreme  case  of  p  =  ±1  the  ellipse  will  degenerate  into  a  straight  line  along 
one  of  the  diagonals  of  the  rectangle.  To  help  visualize  the  concept,  the  bivariate  normal 
distribution  density  function  can  be  thought  of  as  a  surface  above  the  plane.  The  contours 
of  equal  density  are  contours  of  equal  altitude  (equal  probability)  on  a  topographical  map; 
they  indicate  the  shape  of  the  hill  (or  probability  surface)  (Anderson,  1984). 


Figure  2.1  CO VAJUANCE  ELLIPSES 


The  form  that  the  ellipse  of  covariance  was  presented  in  Equation  2.14  did  well  to 
illustrate  the  properties  of  an  error  ellipse.  In  order  to  implement  the  combination 
algorithm  we  must  be  able  to  derive  the  covariance  matrix  from  the  initial  parameters.  This 
is  most  easily  accomplished  by  rewriting  Equation  2.14  in  matrix  form  as 

/=(x-x)''C“’(x-x)  (2.18) 
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where  x  is  a  Gaussian  random  vector,  x  is  the  true  mean  of  the  vector  x  and  C  is  the 
covariance  matrix  of  x  (Roecker,  91).  The  chi-squared  value  with  2  degrees  of  freedom 
corresponds  to  the  probability  that  the  observation  falls  within  the  ellipse  described  by  this 
equation.  In  Equations  2.15  thru  2.17  the  angle  between  the  ellipse  axes  and  Xi,  X2  axes, 
and  the  semi-axes  were  derived  from  Equation  2.14.  The  same  information  can  be 
extracted  from  Equation  2.18  with  the  help  of  the  eigenvalues  and  eigenvectors  of  C 
(Roecker,  1991).  The  eigenvalues  are  found  from 


\e-st\  =  o 

(2.19) 

where  I  is  the  identity  matrix.  Substituting  for  C  from  Equation  2.6  and  taking  the 
determinant  yields 

£^-e{cjf+  a^)  +  (T^al  - =  0 . 

(2.20) 

Therefore,  the  eigenvalues  are 

(t/  +  ±  +  a^y  -  4(ay,  - 

2 

(2.21) 

and  the  semi-axes  are  the  square  roots  of  the  eigenvalues  multiplied  by  the 
freedom  x^ -value  corresponding  to  the  specified  probability  p 

two-degree  of 

(2.22) 

b  =  ylx^S- 

(2.23) 

Values  of  vs.  p  are  given  in  Table  2.1 


7 


Table  2.1  Values  of  Chi-squate  vs.  Ptobabij^ 


f 

P 

2.30 

0.683 

4.61 

0.90 

6.17 

0.954 

9.21 

0.99 

11.9 

0.9923 

18.4 

0.9999 

B.  BAYES’  THEORM 

This  project  is  intended  to  combine  target  error  ellipses  (discussed  in  the  previous 
section)  generated  by  several  independent  observations,  each  of  which  can  be  considered 
Gaussian.  The  result  of  which  will  be  an  error  ellipse  that  represents  the  smallest  area 
containing  the  target  location  with  the  specified  probability  p.  The  nucleus  of  the 
Combination  Algorithm  is  based  on  the  work  of  Nelson  M.  Blachman  of  GTE 
Government  Systems  Corporation.  Since  Blachman’s  application  of  Bayes’  Theorem  is 
fundamental  to  the  solution,  excerpts  of  his  manuscript  (Blachman,  1989)  are  presented 
below; 


We  suppose  that,  on  the  basis  of  I  independent  sets  of  observations  Oj ,  Oj , . . . ,  O;  of 
the  same  target,  a  location  ellipse  has  been  found.  Each  such  ellipse  is  a  contour  of  the 

Gaussian  conditional  probability  density  function  (pdf  )  of  the  target 

location  based  on  one  set  O,  of  observations.  Each  of  the  I  given  ellipses  is  based  on 
the  Bayesian  relationship 
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(2.24) 


Pi(x,y)i 

p(o, 

1 

A 

0) 

expressing  the  probability  density  function  ( pdf)  of  the  target  (x,  j)  conditioned  on  the 
observations  O,  as  the  ratio  of  the  product  of  the  prior  (i.e.  before  O^)  pdf  of 
[x,y)  times  the  pdf  of  0,  conditioned  on  the  location  [x,y)(this  product  being  the 
Joint  pdf  of  and  O^)  to  a  normalizing  constant,  which  is  the  unconditioned  pdf 

of  Oj  and  is  the  integral  of  the  foregoing  product  over  all  x  and  y.  Because  different 
observers  may  utilize  different  prior  pdfs,  Pi{x,y)  can  depend  on  i.  As  long  as 
p.{x,^  is  constant  over  the  entire  area  where  ^0,  it  cancels  out  here 

because  )  is  proportional  to  it,  and  so  its  value  is  not  important. 


For  each  I  the  pdf  of  O,  conditioned  on  target  location 


is  therefore 


(2.25) 


On  account  of  the  statistical  independence  of  the  different  sets  of  observations  for  a 
given  target  location,  the  product  of  Ao\  over  all  i  is  the  conditional  pdf 

of  the  entire  set  of  observations  for  a  given  target  location 
{xyf).  Multiplying  it  by  the  prior  pdf  p{xy^,  \^e  get  the  joint  pdf 
of  the  observations  and  the  target  location.  Dividing  by 
p{p^y027--^0j^,  which  is  the  integral  of  the  latter  over  all  x  and  y,  we  get,  by 
Bayes’  theorem,  the  pdf  of  (x,  conditioned  on  all  observations, 
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(2.26) 


/  i„  A^.y)  ^Ao)A^’)P) 

p(o„...,o;rn 


The  factors  Ao)  and  p{p]^,...,Of  which  depend  only  on  the 
observations,  serve  merely  to  normalize  this  posterior  pdf;  they  do  not  effect  its  shape. 
As  long  as  the  prior  pdf  p{x,^  is  constant  over  the  whole  area  where  the  target  may 

lie,  it  cancels  out  as  before,  and  it  too  does  not  affect  the  shape  of  lfc,y\0„...,0j), 
which  is  thus  determined  as  Equation  2.26  indicates,  by  the  product  of  the  pdfs 
j  based  on  the  individual  sets  of  observations. 


Blachman  has  set  the  theoretical  foundation  from  which  a  practical  combination 
algorithm  can  be  developed.  In  order  to  s^pply  Equation  2.26  to  a  spherical  earth  we  must 
have  the  ability  to  manipulate  the  observations  by  a  constant  scaling  factor.  Roecker  while 
expanding  on  Blachman’s  work  (Roecker,  1991)  noted  that. 


as  long  as  the  prior  pdfs  p{X)  and  p^iX)  are  constant  over  the  areas  of  interest, 
they  will  cancel  out  with  p{0,,,..,0j)  and  p(0,)  to  affect  p[]^0,,...,0j)  only 
as  a  scaling  factor.  This  leaves  the  conditional  pdf  p{x\0^y...^0^  completely 
dependent  on  the  product  of  j  for  its  shape,  and  can  be  expressed  as 

Ax\0,.-.0,)=KYIAAP,)  (2-27) 

i=l 
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where  K  is  the  scaling  factor  and  A'  is  a  vector  containing  the  target  location.  This 
observation  will  later  allow  the  observation  covariance  matrices  to  be  scaled,  compensating 
for  a  spherical  earth. 


C.  METHODOLOGY 

In  order  to  apply  Equation  2.27  a  method  of  extracting  an  optimal  estimate  of  the 
target  location  must  be  developed.  Since  />(x|o,  )  is  Gaussian,  the  conditional  pdf  of  the 
target  location  can  be  written  as 


/?(x|o„...,Oji,)=.S:,nexp|-  (2.28) 

p(il®i,..»M)  =  ^iexp|-^Z(i-o,)^C-'(x-o,)j  (2J29) 


where  Ki  is  the  scaling  factor,  x  the  vector  containing  the  target  location,  o,.  the  i  th 

observation  of  x  and  C  the  covariance  matrix  of  x  (Roecker,  1991).  The  optimal  estimate 
for  any  symmetric  cost  function  is  the  same  for  Gaussian  distributions  (Van  Trees,  1968) 
and  results  by  minimizing  the  pdf  detailed  in  Equation  2.29.,  The  minimum  of  Equation 
2.29  is  found  by  setting  the  first  derivative  to  zero  as  presented  below 


1  dE 


M 


2  ^x 


(2.30) 


where  E^  is  defined  as 


A/ 

1=1 


(2.31) 


and  where 
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(2.32) 


Note  that  M  is  the  number  of  observations  and  n  is  the  number  of  components  of  the 
vectors  x  and  o  (n=2  for  this  problem).  To  begin  with  we  will  examine  the  first  observation 


(2.33) 


hi  ^11  ^ln\ 


-o,„ 


(2.34) 


Z5,i(x,-oJ  -  Z5,„(x,-o,) 


Xi-o„ 


x„  -o, 


(2.35) 


=  Sh; 


(2.36) 


then  by  the  product  rule  we  have 


(2.37) 


since  all  summing  is  to  «  we  can  make  a  change  of  indices  for  readability 


C7X^  ,=1  i=l 


Because  the  covariance  matrix  is  symmetric  we  have 
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(2.39) 


dE, 


n 

1=1 


-o„ 


) 


(2.40) 

In  matrix  notation  this  becomes 

(2.41) 

.  . 

Setting  the  derivative  ^  to  zero  we  have 

2C^’(x-Oi)  =  0 

(2.42) 

ft 

1 

ft 

li  . 
o 

(2.43) 

X  =  0^ 

(2.44) 

So  the  initial  estimate  for  x  is  the  first  observation  o, ,  which  was  to  be  expected 
Now  that  the  process  for  a  single  observation  has  been  established,  we  can  generalize  to 
multiple  observations.  To  begin  with  we  expand  the  reduction  variable  E  to  include  M 
observations 


(2.45) 


and  once  again  set  the  derivative  to  zero 
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m~\ 

M 

m=l 

f  M  '\  M 

Z«;' *-E(«>.)='>  p'«) 

Vm=l  /  m=l 

M 

c-i-l€(e>.)=o  (2.50) 

m-l 

Which,  with  a  little  algebraic  manipulation  yields 

M  f  M 

x  =  Cf](C>„)  where  *=2®;’  (2.51) 

m=l  ^««=1  ^ 


As  depicted  in  Figure  2.2,  x  is  the  vector  representing  the  location  of  the  weighted 
average  of  M  number  of  o  observations  and  C  is  the  resultant  covariance  matrix.  This  is 
die  optimal  estimate  of  the  target  location  we  desire  and  the  premise  upon  which  the 
Combination  Algorithm  is  based. 


Figure  2.2  OPTIMAL  ESTIMATE  OF  TARGET  LOCATION 
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III.  CONCEPT  DEVELOPMENT 


A.  COVARIANCE  MATRICES 

1,  Scaled  Covariance  Matrix 

Chapter  II  laid  the  theoretical  groundwork  for  the  Combination  Algorithm.  We 
have  seen  that  the  target  error  ellipses  we  are  attempting  to  combine  can  be  modeled  by  an 
ellipse  of  covariance.  The  covariance  matrix  being  the  structure  defining  the  target  error 
ellipse.  In  the  real  world  problem,  the  error  ellipse  will  have  been  generated  upon  the 
spherical  earth.  In  order  to  reduce  the  overall  problem  to  a  two  dimensional  one,  each 
ellipse  will  be  projected  into  a  plane.  Therefore,  the  covariance  matrices  will  need  to  be 
modified,  compensating  for  this  projection  from  the  spherical  earth  to  the  plane.  This 
modification  is  to  compensate  for  the  inequality  between  latitude  error  and  longitude  error 
as  the  latitude  increases  from  zero.  The  geometry  of  the  projection  from  a  spherical  earth 
to  a  tangent  plane  is  depicted  in  Figure  3.1.  We  need  to  find  the  distance  (length  of 

projected  ellipse  axes)  when  given  the  arc-length  along  the  surface  of  the  sphere  ci,pi,ere 

ellipse  axes),  where  =  radius  of  earth  and 

aspHere=red  (^-l) 

0  =  (3.2) 


Observing  the  geometry 


tan»  =  ^  (3-3) 

apiane=r,imd  (3.4) 


which  leads  to  the  transformation 
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Figure  3.1  PROJECTION  GEOMETRY 


The  congruence  transformation 

C  =  S^CS  (3.6) 


produces  the  scaled  covariance  matrix  C  that  has  maintained  the  symmetry  of  the  original 
covariance  matrix..  The  scaling  matrix  S  is  obtained  by  observing  from  Figure  3.2  below 
that  the  arc-length  measured  along  the  surface  of  the  earth  at  constant  latitude  (f) 

corresponding  to  a  change  in  longitude  is  given  by  cos(^)AA  ;  and  the  arc-length 
measured  along  the  surface  of  the  earth  at  constant  longitude  A  corresponding  to  a  change 
in  latitude  A^  is  given  by  To  put  latitude  and  longitude  on  equal  footing  (in  terms 

of  arc-length)  changes  in  latitude  are  multiplied  by  and  changes  in  longitude  are 

multiplied  by  cos(^i)  .  Therefore, 


(3.7) 
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A..  1 

K/ 

cP  X  \ 

\ 

/  AL  \ 

\ 

/  \ 

y\^  1 

Figure  3.2  GEOMETRY  OF  SCALING  MATRIX 


Expanding  Equation  3.6  we  find  that 


c 


0  ^ 
vO  r^cos<j>j 


r 


^f¥ 

°  1 

lo 

cos^j 

(r. 

0  ^ 

COS^j 

u 

cos^j 

f  1  2  j  ^ 

^  cos^Zi 


(3.8) 


(3.9) 


(3.10) 


This  is  the  scaled  covariance  matrix,  where  (j)  =  latitude,  X  =  longitude  and  r^  is  the  radius  of 
the  earth.  In  the  section  to  follow  it  will  be  necessary  to  solve  for  the  eigenvectors  of  this 
matrix.  To  aid  in  that  process,  the  following  notation  will  be  used  to  refer  to  the  scaled 
covariance  matrix  in  Equation  3.10. 


P) 


(3.11) 
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where  use  has  been  made  of  the  fact  that  a ^ 


2.  2x2  to  3x3  Covariance  Matrix  Conversion 

The  objective  is  to  convert  the  2x2  covariance  matrix  given  in 
longitude  (Equation  3.11)  to  a  3x3  covariance  matrix  in  terms  of  (x,y,z) 
is  as  follows 


where  the  3x3  covariance  matrices  are 


r  __  ^ 


O^xz 

c,  ,  = 

{x,y.2) 

^yy 

^zz) 

^(jiX 

\ 

^hX 

^hhJ 

and  the  3x3  transformation  matrix  T  is  given  by 


'  6x 

dx 

dy^ 

d(t> 

dh 

dy 

dy 

d(t> 

dX 

dh 

dz 

dz 

dz 

~dhj 

where  for  a  spherical  earth 

X  =  +/j)cos<^  cos  A 

><  =  (r^ +/j)cos^  sin  A 


terms  of  latitude  and 
.  The  transformation 

(3.12) 

(3.13) 

(3.14) 

(3.15) 

(3.16) 

(3.17) 
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(r^  +h)sin^ 


Te  being  the  radius  of  the  earth  and  the  partial  derivatives  are 


dx 


-  )cos(/l ) 


dx 

Ja 


~(re  cos(^  )sin(/l ) 


—  =  cos(^)cos(A) 


d(l> 


=  -(r^  +  h)  sin(^  )sin(A ) 


Jj=(r,  +/»)cos(^)cos(2) 


dy 

—  =  cos(^)sin(A) 
ah 


dz 

'Ji 


=  +^)  cos(^J) 


dz 

~dX 


=  0 


(3.18) 

(3.19) 

(3.20) 

(3.21) 

(3.22) 

(3.23) 

(3.24) 

(3.25) 

(3.26) 

(3.27) 


To  define  we  start  by  finding  the  unsealed  covariance  matrix  from 


Equation  3.6 


(3.28) 
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where 


S"' 


n 

0 


0 

1 

r,cos(^j 


(3.29) 


and  the  scaled  covariance  matrix  C  is  given  in  Equation  3.11.  Using  the  unsealed 
covariance  matrix  in  terms  of  latitude  and  longitude  and  assuming  the  hei^t  h  to  be 


zero  with  zero  uncertainty,  we  can  form  the  matrix 


K  0 


0  Oj 


(3.30) 


3.  Selection  of  Eigenvalues  and  Eigenvectors 

The  key  factors  in  the  transformation  of  an  observation  seen  in  terms  of  semi-axes 
a,  b  and  heading  0,  to  it’s  description  by  a  covariance  matrix,  are  it’s  eigenvalues  and 
eigenvectors.  This  section  will  accomplish  four  tasks;  derive  the  eigenvalues  from  the 
characteristic  polynomial  of  Equation  3.11,  derive  various  identities  from  the  characteristic 
polynomial  that  will  be  useful  in  later  sections,  normalize  and  solve  for  the  eigenvectors,  and 
determine  which  of  the  solutions  will  be  numerically  stable. 

a.  Eigenvalues 

In  order  to  produce  the  eigenvalues,  which  will  be  denoted  by  s,  the 
nonlinear  equation 


(C-£l)x  =  0  (3.31) 

needs  to  be  solved.  To  begin  the  process  the  scaled  covariance  matrix  C  is  shifted  by  ffl, 
where  I  signifies  the  identity  matrix 
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C-£l  = 


(3.32) 


a-e  Y 
.7 


By  definition  (Strang,  1988)  the  number  e  is  an  eigenvalue  of  C  if  and  only  if 
det(C  -  fl)  =  0 .  This  determinant  leads  to  the  characteristic  polynomial.  It’s  roots,  i.e.  where 
the  determinant  is  zero,  are  the  eigenvalues. 


det(C-£l)  =  (a-5)(y0-e)-r' 

=  s^-(a  +  p)e  +  ap-r^  (3-34) 

The  real  solutions  of  the  characteristic  polynomial  (the  eigenvalues)  are  found  by  the 
quadratic  formula  as  show  below. 


= 


(a  +  /3}±/{  a  +  fiY  -4afi  +  4Y^ 
_  _ 


(3.35) 


(g  +  (3) ± +  2ap  +  - 4ap  +  4y'^ 

_  - 


(3.36) 


(a  +  yS) -2aP  +  +4^’ 

2 


(3.37) 


(a  +  P)±V(a-P)^+4Y^ 


(3.38) 


b.  Identities 

With  the  characteristic  polynomial  at  hand,  this  is  a  good  time  to  observe 
two  identities  that  will  become  useful  in  future  sections.  Both  of  the  following  identities  are 
obtained  from  factorization  of  the  general  quadratic  equation,  which  is 
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[s-s,\s-s_)  =  0 


(3.39) 


-[s^  +s.)e  +  s^s_  =0. 


(3.40) 


Comparing  Equation  3.40  with  the  characteristic  polynomial  Equation  3.34  we  find  two 
identities  that  will  later  aid  in  establishing  relationships  between  the  scaled  covariance  matrix 
Equation  3.11  and  the  primary  data  parameters  a,  b,  and  0, 


s^+s_  =a+P 


(3.41) 


s^s_  = 


Also  worth  noting,  but  of  less  significant  impact,  are  the  identities 


(3.42) 


(3.43) 


=  a  and  s_  .=  P  ■ 

+  r=o  r=o 


(3.44) 


= 


c.  Eigenvectors 

This  section  will  present  solutions  for  the  eigenvectors  x,  where 
,/  =  ± .  Writing  Equation  3.31  in  terms  of  Equation  3.32  and  the  notation  just 


mentioned  yields 


^a-6  y  Y^i^ 


\  y  p-e) 


(3.45) 


Equation  3.45  represents  a  system  of  two  equations  with  two  unknowns.  By  adding  the 
formula 


xl  +XI  =1 


(3.46) 
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to  the  system  of  equations,  we  can  simultaneously  normalizing  the  eigenvector  to  unit 
length  and  supply  an  additional  formula  that  will  be  used  to  produce  four  valid  normalized 
solutions.  Each  pair  of  equations  will  generate  solutions  for  the  positive  and  negative  roots 
of  Equation  3.38.  Each  solution  will  be  checked  for  computational  stability.  The  most 
stable  vector  from  each  root  will  eventually  be  used  to  populate  the  eigenvector  matrix. 

(1)  Solution  for  .  The  first  solution  presented  is  for  x+^ , 
which  corresponds  to  the  solution  of  the  first  equation  in  Equation  3.45.  When 

{a  -  )x,  +yx2=0  (3.47) 


x,= - (3.48) 

a-s^ 

The  second  equation  needed  to  solve  this  system  is  supplied  by  normalizing  the 
eigenvectors  to  unit  length  using  Equation  3.46.  By  substitution  we  have 


Solving  for  we  get 


Substituting  Equation  3.51  back  into  Equation  3.46  yields 


-r 


^j(a-s^y  +r' 


(3.49) 


(3.50) 


(3.51) 


(3.52) 
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Therefore, 


Check: 


(1)  ^ 


2  2 

+r 


r  _y  ^ 

'va  -  s^J 


1 

1 

pAa-sJ 


^ -  ay  +  aY  -  ys^'^ 


\£^s_  -PsJ 


r  \ 


+  r 


Now  Equation  3.41  states 


s_- P  =  a-s^  , 


therefore, 


and  Equation  3.31  holds  true. 

In  checking  for  numerical  stability  note  that  when  y-^O  we  obtain 
of  x+^ ,  therefore  is  not  numerically  stable  as  y->0. 


(3.53) 


(3.54) 


(3.55) 


(3.56) 


(3.57) 


for  the  components 


24 


(2) 


(2)  Solution  for  xf  ^ .  The  second  solution  evaluated  is  x^  , 
which  corresponds  to  the  solution  of  the  second  equation  in  Equation3.45  when  S-S^ 


yx,+[p-s^)x^  =0 


X,  = 


r 


(3.58) 


(3.59) 


The  second  equation  needed  to  solve  this  system,  as  in  the  previous  solution,  is  supplied  by 
normalizing  the  eigenvectors  to  unit  length  using  Equation  3.46.  By  substitution  we  see 


xj  =  1 


(3.60) 


solving  for  Xj  we  get 


x^  = 

■*•2 


^2 


^{p-sYY 


(3.61) 


(3.62) 


Substituting  Equation  3.62  back  into  Equation  3.46  yields 


X,  = 


s,-P 


+y' 


(3.63) 


Therefore 


= 


1 

\  Y  ) 

(3.64) 
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Check: 


Cx 


(2)  _ 

+ 


+r 


y^ 

\y 

PJ 

\  y  ) 

P- 


as^  -aP  +  y 


2^ 


(3.65) 


(3.66) 


-s^8^ 


}I(P~^^) 


2  2  V  ys^  J 


+y 


V(^ 


+y 


^  a  - 
\  y  ) 


Using  the  Equation  3.41 


a-s_  =  e^-  P 


yields 


Cxl^>=£X^> 

and  Equation  3.31  holds  true. 

In  checking  for  numerical  stability  we  note  that  when  y— >0 


1 

y 

r-(|3-a)) 

=  f'l 

a-(3 

0  J 

UJ 

(3.67) 


(3.68) 


(3.69) 


Hence,  is  a  numerically  stable  solution.  Therefore  we  will  use  Equation  3.64  for  the 
solution  from  the  eigenvalue  .  The  next  step  is  to  determine  which  of  the  two  equations 
from  the  eigenvalue  S_ .  will  be  used. 
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(3)  Solution  for  .  The  next  solution  presented  is  for 

which  corresponds  to  the  solution  of  the  first  equation  of  Equation  3.45  with  £  = 

(a-£_)x,  +yx2=0 

3 

(3.70) 

II 

1 

I 

(3.71) 

The  second  equation  needed  to  solve  this  system  of  equations  is  supplied  by  normalizing 
the  eigenvectors  to  unit  length  using  Equation  3.46.  By  substitution  we  see 

+  1  X,  -  1 

J 

(3.72) 

Solving  for  x^  we  get 

^2  (  V  ,  2 

\a-e_)  +/ 

(3.73) 

(a-e_) 

^l(a-e_f  +r^ 

(3.74) 

Substituting  Equation  3.74  back  into  Equation  3.46  yields 

yl{a-e_y+r^ 

(3.75) 

Therefore, 

1  f  ) 

(1) 

li  V  ,  .,2  Voc-eJ 

}l{oc-s_)  +r 

(3.76) 
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Check: 


Cx 


(1)  _ 


^l(oc-£^) 

2  2I 
+r 

1 

( 

)■ 

^  2I 

+  ^2  V. 

1 

/ 

yl(a-£.) 

2  2  1 
+r 

£_ 

i 

^  a  y^ 

'  -r  ^ 

vr  p) 

\a  -  £_J 

ay  +  ay  -  ys_ 

.,2  ,  ^  R„ 


-ys_ 
\^s^s_  -fisj 


+  y 


-r 


(377) 


(3.78) 


(3.79) 


(3.80) 


Now  Equation  3.41  states 


£^-  P  =  a-s_ 


Therefore, 


showing  that  Equation  3.31  holds  true. 

In  checking  for  numerical  stability  we  note  that  when  y->0 


1 


a  -  P  va  -  pj 


fo\ 


(3.81) 


producing  a  stable  solution. 
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(2) 

(4)  Solution  for  .  The  last  solution  presented  is  for  xi  , 
which  corresponds  to  the  solution  of  the  second  equation  in  Equation  3.45  when  s- 


[p-e.)x^=0 


(3.82) 


X,  = 


P-S- 

r 


(3.83) 


The  second  equation  needed  to  solve  this  system,  as  in  the  previous  solutions,  is  supplied  by 
normalizing  the  eigenvectors  to  unit  length  using  Equation  3.46.  By  substitution  we  see 


[Izll 


+1 


1X,^=1 


(3.84) 


Solving  for  we  get 


7 


■\-Y 


(3.85) 


(3.86) 


Substituting  Equation  3.86  back  into  Equation  3.46  yields 


^(p-s.f  +r" 


(3.87) 


Therefore, 


x®  = 


Ap-^-T 


+  7 


'^s_  -  P^ 
\  7  7 


(3.88) 
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Check: 


Cx 


(2)  _ 


+  r‘ 


a  Y 

r  P) 


fi 


V  r 


T\ 


f  as _  -aP  +  Y 

^  -Py^Py^ 


^(p-^-f 


+y 


^ as _  -  s^s^ 

YS- 


{3m) 


(3.90) 


(3.91) 


(3.92) 


Using  Equation  3.41 


a-  s^=  s_-  P 


in  Equation  3.92  gives 


showing  that  Equation  3.31  holds  true. 

Checking  for  numerical  stability  we  note  that  when  y->0  we  obtain  ^  for  the  components 

of  ,  therefore,  is  not  numerically  stable  as  y->0. 

The  eigenvalues,  like  the  eigenvalues^,  produced  one  computationally  stable 
eigenvector.  Note  that  s_  >  0  ands^  >  0  so  both  values  are  positive.  The  equations  to  be 
used  in  solving  for  the  eigenvector  matrix  G  wiU  be:  Equation  3.64  which  corresponded  to 


30 


the  solution  of  the  second  equation  in  Equation3.45  when  S  =  8^,  and  Equation  3.76  which 
corresponds  to  die  solution  of  the  first  equation  of  Equation  3.45  with  s-s_. 


d.  Eigenvector  Matrix 

This  section  will  develop  the  solution  for  eigenvector  matrix  G  and  verify  it 
by  producing  die  diagonalized  eigenvalue  matrix  A.  The  eigenvector  matrix  G  is  the  matrix 

whose  columns  are  the  eigenvectors  and  (Recall  these  are  die  computationally 
stable  eigenvector  solutions  found  in  sub-section  c.  of  this  chapter).  In  order  to  simplify 
the  construction  of  G,  some  algebraic  manipulation  of  the  vectors  is  required.  Consider 
Equation  3.64.  Since  s^-  P  =  a-8_  (Identity  Equation  3.41)  we  have 


j(2)  ^ 


1  1 

^[a-s_\s^  -P)  +  A 

\  r  ) 

(3.93) 


Now  consider  Equation  3.76.  Again,  since  a  -  £_  =  fi;  -  we  can  say  that 


1 

(  -y  \ 

^{a-s_\s^-p)  +  A 

la  -  ej 

By  combining  Equation  3.93  and  Equation  3.94  die  eigenvector  matrix  G  is 

'  p)^  A  '' 


(3.94) 


(3.95) 


The  determinant  (difference  of  the  products  of  the  numbers  in  die  two  diagonals)  of  G  is 


det(G)  = 


^{a-e.\e^-p)  +  y‘ 


\e,-p'ia-s.)  +  r^},  (3.96) 


which  is  equal  to  one,  and  the  inverse  of  G  is 
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G’  = 


^a-s_ 


yl{a-£.X^^-p)  +  r' 


■r 


(3.97) 


The  diagonalized  eigenvalue  matrix  A(the  diagonal  matrix  with  the  eigenvalues  and 
the  diagonal  components)  is 


as 


G  ’CG 


[a-s_)[e^-  p)  +  r 


s^-P  -r 


a  Y 

\-r  P) 


a -s_  Y  I 

^  -Y  e^-PJ 


k  -p] 

1 

{a-s.) 

{e,-p)  +  Y^ 

1 

{a-s_] 

(s.-zsl  +  r' 

1  1 

Y  a-eJ^ 

a-s_  Y  yias^-ap^Y^  -UY+aY-Y^. 


a-s_ 


a-s_  Y 


ae^-s^e_  -Y^- 


'eXa-s_)  -Ye-  ^ 
Ye,  eXe.-p)} 


(3.98) 


(3.99) 


(3.100) 


(3.101) 


s,{a-e_)  +Y^e^  -  Ye-{ec  -  e.)  +  Ye-{e+  p) 


(3102) 


and  by  using  a-s_  =  s^-p  (Identity  Equation  3.41)  this  becomes 


A  = 


(e,  0^ 
vO  e_j 


(3103) 


With  this  check  it  is  reasonable  to  assume  that  Equation  3.95  is  in  fact  the  eigenvector 
matrix. 


4.  Variable  Relationships 

In  die  preceding  text  there  have  been  two  sets  of  relationships  established.  In 
Selection  of  Eigenvalues  and  Eigenvectors  (Chapter  III.A.2)  a  relationship  between  the 
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elements  a,  P  and  y  of  the  scaled  covariance  matrix  Cand  the  eigenvalues  5^.  ,  s_  was 
established.  A  second  relationship  between  the  eigenvalues  and  the  data  parameters  a,  b,  0 
and  was  detailed  in  Error  Ellipse  Relationship  to  the  Bivariate  Normal  Distribution 

(Chapter  II.A).  The  following  section  uses  these  relationships  to  develop  C  in  terms  of  a, 
b,  0  and 

The  eigenvector  from  the  positive  root  contains  information  relative  to  the 
direction  of  the  semi  major  axis  of  the  ellipse.  Of  the  two  positive  root  solutions,  it  was 
determined  that  the  second  was  the  more  numerically  stable.  Therefore,  the  x,  y 

components  of  the  eigenvector  can  be  used  to  determine  0  as  illustrated  in  Figure  3.3. 


From  Equation  3.64  comes  the  following 


tan^  = 


X 


(2) 

2+ 

.(2) 

M+ 


(3.104) 


The  semi-axes  for  the  ellipse  corresponding  to  a  probability  level  p  are  related  to  the 
eigenvalues  and  s_ .  By  Equations  2.22  and  2.23  we  see  that 
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(3.105) 


and 


where  x"  is  the  two-degree  of  freedom  value  of  corresponding  to  probability  level  p. 
Manipulating  Equation  3.104  we  find 

-  P  =  y  cot0  (3.107) 


which  leads  to  the  formula 


P  =  e^  -ycoXO . 

(3.108) 

From  the  Identity  Equation  3.41  it  follows  that 

a  =  s^+s_-P 

(3.109) 

=  e^+e_-s^  +  ycoXd 

(3.110) 

=  s_  -^-ycoXO. 

(3.111) 

Substituting  Equation  3.108  and  Equation  3.111  into  the  Identity  Equation 
second  order  polynomial 

3.42  we  find  the 

s^s.  =  (^.  -;^cot^)-r' 

(3.112) 

=  —  ye^  cot^  +  ye^  cotd—y^  cot  0  —  y 

(3.113) 

This  can  be  rewritten  as 
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-  (cot^  0 + 1)  +  r(5+  cot Q-S-  cot  ^)  =  0 

(3.114) 

or 

Y ^ (cot^  ^  + 1)  +  y{^-  cot 0-s^cot0)  =  O. 

(3.115) 

Y {^'[cot^  ^  + 1] + [e,  cot  0-s^  cot  ^]j  =  0 

(3.116) 

The  non-2ero  solution  of  this  equation  is 

cot  0-s_  cot  0 
co\}0  +  \ 

(3.117) 

COt^ff+  -f_) 

1 

(3.118) 

sin^  0 

=  sin^  ^cot^fi;  -  s_). 

(3.119) 

Seen  in  a  more  convenient  form  throu^  application  of  the  trigonometry 

identities  as 

Y  =  sin  ^cos  0{e^  ~  • 

(3.120) 

Recalling  Equations  3.105  looking  for.  Back  and  3.106  we  see  that  Equation  3.120  is  y  m 

terms  of  a,  b,  0  and  yC  >  exactly  what  we  were  substituting  for  a  and  p  we 

:  have 

a-6_  +^cot0 

(3.121) 

=  g_  +  sin^cos0(e;  -  e^coX0 

(3.122) 

=  8_+  COS^  0S^  -  COS^  06_ 

(3.123) 

leading  to  the  useful  fonn 
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a  =  e_+  cos^  d{s^  -  e.) . 

(3.124) 

For  P  we  get 

(3.125) 

=  s^-  sin0cos^e+  -  e_)cat9 

(3.126) 

leading  to  the  form  we  want 

P=  8^-  cos^  -  e_). 

(3.127) 

Equations  3.120,  3.124  and  3.127  establish  a  direct  relationship 
covariance  matrix  and  the  data  parameters  a,  b,  0  and  %  . 

between  the 

B.  N+1  COMBINATORIAL  EFFECTS 

The  ability  to  combine  multiple  error  ellipses,  yielding  the  position  vector  and  the 
associated  covariance  matrix  is  provided  throu^  Equation  2.51.  This  section  is  concerned 
with  the  effect  of  including  an  additional  error  eUipse  to  the  existing  solution,  i.e.  the  N+1 
case.  In  order  to  more  clearly  describe  this  problem,  let  us  look  at  an  example  containing 
three  ellipses.  The  question  posed  is  will  the  solution  for  the  case  when  all  three 
observations  are  taken  at  once,  as  depicted  in  Figure  3.4, 
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be  equal  to  the  solution  for  the  case  in  Figure  3.5,  where  two  elUpsoids  are  initiaUy 
processed,  then  a  third  is  combined  to  that  solution? 


r\ 

A 

u 

cr-  ^ 

11 

^  1 

I 

1 

Figure  3.5  COMBINATION  EXAMPLE  (b) 


Formally  stated  in  Equation  3.128,  we  need  to  prove  that  the  probability  of  the 
solution  vector  x  given  N+1  observations  is  equivalent  to  the  probability  of  x  given  N 
observations  given  an  additional  observation. 


p(p(xlo,...o^)|o^^,)  =  P(x|o,...o^„)  (3.128) 

The  optimal  estimate  for  the  combination  of  N  observed  error  ellipses  was  derived  in 
Chapter  II.  Substituting  the  scaled  covariance  matrix  Cinto  Equation  2.51  we  see  that 
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X  = 


(3.129) 


.  -1  > 


Se-M  =«;  !«.■'» 

Vf=i  j  V,=1  J  \i=\ 


where 


N  A 


-1 


«;=  S«r'o,  ^  er  =  S«r‘ 

v,=i  y  i=\ 


(3.130) 


C' is  3.  new  variable  btroduced  to  differentiate  between  the  previously  summed  and  the  yet 
to  be  summed  N+1  cases.  The  scaled  covariance  matrix  for  the  solution  of  N+1 
observations  is 


cr.', +ei'J =(cr'  +cj.,) 


(3.131) 


(3.132) 


Now  let’s  work  througji  the  N+1  solution 


(3.133) 


tc:*+C-0  ficro,+C-*,o^J  (3.134) 
y  v,=i  y 


(3.135) 


fe-'  +C-V,)x^,.  =  ic:‘o,  +C-V,o^,, 

'  1=1 


(3.136) 


where  Clf  =  I 


'  ■'■®ArVi)*Ar+t  “  *JV  ■*‘®w+l®Ar+l 


(3.138) 


‘^AT+l 


=(cr  +C;'., )■'[«;-'!„  +e;'.,o,.,]  (3,139) 


EC)>o, 


1=1 


(3.140) 


Which  agrees  with  Equation  3.128  proving  the  N+1  case  is  valid.  This  proof  is  critical 
from  a  software  engineering  point  of  view.  Consider  the  N+1  case  not  valid.  A  piece  of 
software  wanting  to  keep  a  track  updated  over  some  amount  of  time,  would  need  to  store 
the  information  from  all  prior  observations  in  order  to  update  it’s  location  and  associated 
error  ellipse.  Since  the  N+1  case  is  valid,  only  the  most  recent  solution  needs  to  be 
retained.  Allowing  all  previous  fixes  to  be  discarded,  freeing  memory  for  other  tasking. 


C.  ELUPSE  ROTATION 

A  change  of  axes  will  be  required  in  order  to  place  the  error  ellipses  that  are  to  be 
combined  in  an  approximately  orthogonal  coordinate  system.  In  the  latitude/longitude 
coordinate  system,  as  latitude  increases  the  coordinate  system  becomes  progressively 
distorted.  The  closest  approximation  to  an  orthogonal  coordinate  system  within  the 
latitude/longitude  coordinate  system  is  found  by  transforming  the  ellipse  center  to  a 
coordinate  system  where  X,=0.  Where  <j>  is  latitude  and  X.  longitude.  This  coordinate 
system  will  be  referred  to  as  the  prime  or  coordinate  system. 

The  process  will  begin  by  defining  the  coordinate  system  axes  in  terms  of 

^  and  X  with  respect  to  the  original  ex,ey,e2  coordinate  system.  Figure  3.6  depicts  all 
axes  and  angles  relevant  to  the  transformation  between  the  two  coordinate  systems. 
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Figure  3.6  ROTATION  COORDINATE  SYSTEM 


With  die  aid  of  Figure  3.6  the  axis  of  die  prime  coordinate  sptem  e^,  and 
are  defined  as  follows 

=  cos(X)cos((l))e;c  +  sin(X)  cos(<l>)e^  +  sin((|))e^  (3-141) 

=  cos(it  +  X)  ~  +  sin(X  +  tc)  (3.142) 

=  -  cos(A)  sin(<^)e,  -  sin(A)  sin(^)e^  +  cos(^)e^  (3.143) 

As  a  simple  check  we  verify  that  and  are  orthogonal  by  confirming  their  dot  product 
is  zero. 


=-cos^(X)sin(<|>)cos((l))-sin^(X)sm(<l>)cos((l>)  +  sin((l>)cos((l))  (3.144) 


The  third  axis  is  easily  found  by  the  cross-product  of  the  two  previous  directional  vectors 


=  -sin(X)cos(X)sin((|))cos(<l>)e2  +cos(X.)sin^((t))e^  -i-sin(^)cos(X,)sin((l»)cos(<|))e2 


-  sin(X)  sin^  ((|))ej^  +  cos(X)  cos^  (<|))«^  -  sin(X,)  cos^  (3.147) 

=  -  sin(X)e;e  +  cos(A,)e^  (3.148) 

To  ensure  thorou^ness,  the  remaining  combinations  of  axes  are  checked  to  be  orthogonal 
e^-e^  =-  sin(X)cos(X)cos((j))  -t-  sin(X)cos(X)cos((|))  =  0  (3.149) 

•  eri  =  sin(X)  cos(X)  sin(<|))  -  sin(X,)  cos(X)  sin(<t))  =  0  (3.1 50) 


The  unit  vectors  in  the  coordinate  system  can  be  express  as  a  relationship 

of  the  directional  cosines  relative  to  the  e^,ej^,eQ  coordinate  system.  In  unit  terms 


—  cosilix  >  which  leads  to  the  following 

•  %)%  (3-^51) 

=  cos(X.)cos((j))e^  -  sin(A,)«Ti  -  cos(X)sin((|))e^  (3.152) 

=  (^y  ■  h)h  ^[^y ' +  [^y ' 

=  sin(X)cos(<l))e^  +  cos(X.)eTi  -  sin(A,)sin(<l))e^  (3.154) 

4  =  (4  •  +  (^z  -^11)^11+  (4  •  %)k 
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(3.156) 


=  sin(<|>)e^  +  cos((|))e(; 


If  we  have  coordinates  x  = 


y 

\zj 


in  the  coordinate  system,  the  coordinates  of  the 

same  point  in  the  ^incoordinate  system  are 


f  A 

^  A. 

fj\ 

e.  •  0 

X 

1  = 

n 

= 

A.  A 

A  A 

e  •  e 

V  y 

A  /V 

e  *e 

y 

C  y 

\z) 

^  cos(X.)cos((|)) 

sin(X.)cos(<t>) 

sin(<l>)^ 

fx^ 

-  sin(X) 

0 

0 

0 

y 

^-cos(A,)sin(<|)) 

-sin(X)sin(<t>) 

cos(<l>); 

\z) 

The  inverse  transform  is 


A 

e  •  e. 

A  A 

e  -e 

A  A 

A  A 

e  • 

y 

U; 

e  *  e 
y  7 

e  •  e 

e  -Sr 
y  C 

-v 

n 

Id 

^cos(X,)cos((j>)  -sin(A,) 
sin(A,)cos((l>)  cos(X) 
<  sin((l))  0 


-cos(X,)sin(4))^ 
-  sin(X,)  sin((j)) 
cos((l))  j 


'^1 

Id 


(3.157) 


(3.158) 


(3.159) 


(3.160) 


If  we  let  (|)i  ,  Xi  denote  the  latitude  and  longitude  of  a  point  in  the  coordinate 

system.  Then  x^  can  easily  be  expressed  in  terms  (|»i  and  A.i 


= 

Ti 

^cos(^) 

cos(<^,) 

COs(>^5)^ 

lsin(A,) 

U,j 

^  sin 

U)  J 

(3.161) 
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The  coordinates  of  the  point  in  the  coordinate  system  are 


^  cos(A)cos(^)  sin(/l)cos(^)  sin(^P  cos(^  )cos(>^,)^ 

^,  =  =  -sin(A)  cos(i)  0  cos(^.)sin(A,)  (3.162) 

V-cos(A)sin(^)  -sin(A)sin(^)  cos(^)jt  sin(^,)  ; 

If  denote  the  latitude  and  longitude  of  this  point  in  the  coordinate 

system,  then 

^cos(^)cos(/l5)^ 

^.  =  7;  =  cos(^)sin(iJ^)  .  (3.163) 

UJ  t  sin(^)  ; 

Equation  3.163  can  be  used  to  derive  the  latitude  and  longitude  of  the  point  in  the  prime 
coordinate  system.  Table  3.1  shows  the  range  or  value  of  the  longitude  of  the  point  after  it 
is  rotated. 


$,  =sin^(C/)  (3.164)  < 
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In  order  to  transpose  the  heading  of  an  error  ellipse  between  die  two  coordinate 
systems  we  must  know  the  relationship  between  the  semi-major  direction  vector  and  the 
north  pole.  The  coordinates  of  the  north  pole  in  the  Cx,&y,^z  coordinate  system  are 

The  coordinates  of  this  point  in  the  coordinate  system  are  found  when  Equation 

3.158  is  applied  to  the  north  pole  from  the  coordinate  system 


np 


0 

0 

1 1 ) 


Inp- 


f 

Vnp 

■=. 

\CnpJ 

\ 

cos(A)  cos(^) 
-  sin(A) 
-cos(A)sin(^) 


sin(A)  cos(<^) 
cos(A) 

-  sin(A)  sin(^) 


sin(^P 

^0^ 

0 

0 

cos(^); 

a. 

(3.167) 


0 

Uos((l))J 


(3.168) 


Using  Equation  3.168  we  find  the  latitude  and  longitude  of  the  north  pole  in  the  prime 
coordinate  system 


^np  =  sin  ^(cos((l))) 

(3.169) 

Cin(4|)] 

r0if(|>>0 

[71  if  (|)<  0 

(3.170) 

D.  ELUPSE  HEADING  RELATIVE  TO  NORTH  POLE 

The  individual  error  ellipse  headings  relative  to  true  north  are  changed  during  the 
axes  transformation.  Consider  an  eUipse  at  some  arbitrary  latitude,  with  a  heading  of  true 
north.  When  the  ellipse  is  transformed  to  the  prime  coordinate  system,  the  heading  is  no 
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longer  valid  as  depicted  in  Figure  3.7.  This  section  will  present  the  process  which  calculates 
this  change  in  heading. 


The  first  step  in  this  process  is  to  determining  each  error  ellipse  heading  relative  to 
true  north  .  Let  two  points  have  latitude  and  longitude  and  We  want  to  find 

the  heading  from  A  to  B  at  point  A  measured  relative  to  true  north.  The  following  pseudo 
code  addresses  the  four  trivial  cases:  one  of  the  points  is  located  on  either  of  the 
poles,  =  Xg  and  X^=^  X^±  tv  .  First  a  point  located  on  the  north  pole 


^  *I>B  then  (o^=TV 

If  ^  f  then  =  0 

For  the  south  pole 

If^  =-f»  4  ^ then  =  0 

If  =  -|>  -f  then  Oab  =  ^ 

Now  the  case 


If  4  =  4 

^  A 

<0^=0 
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Else  If 

and  the  final  trivial  case,  =  Xg±7t 

IfA,  =4±;r 
If  >  -<f>B 
=0 

Else  If  <-<l>B 

G)ab=^ 

If  no  trivial  cases  are  valid  apply  Napier’s  analogies. 

With  the  trivial  cases  out  of  the  way  we  move  on  to  examine  the  cases  where 
Napier’s  analogies  will  need  to  be  applied.  Refer  to  Figure  3.8  for  notation  definitions  used 
in  die  following  problems.  First  consider  the  case  where  the  longitude  of  point  B  is  greater 
than  that  of  point  A 


AX  =  Xg  — 

AX  >  0  ;  AX  <  7C 

71  ,  ,71 

<1)^  ;t±-  ;  (|)5  ?^±- 


N 

Va  I 

C'X'^AB 

^  /  1 

s 

B 

Figute3.8  HEADING  BETWEEN  PODJTSBANDA 


The  points  bvolved  in  this  problem  Ue  upon  the  surface  of  the  earth.  Taking  this  into 
account  we  must  look  to  spherical  trigonometry  for  a  solution.  Napier’s  analogies  (Selby, 
1967)  provide  the  relationships  from  which  this  problem  can  be  solved.  Employed  below 

we  find  that 
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(3.171) 


[o  ^5  -  ©  BA  ]]  -  cot(^  2  f  1 


^  sin 

[i 

--r  +  ^A 

l2 

_2  2  J/ 

sin 

f\ 

U 

sin[^[<|)^-(t>5] 

cotl-AX' 

\2 


sinQ[7C-<t)^  -^B]j 


(3.172) 


\  sinl 

cotl  —  AX,) — T 
<.2  /  [ 
coa  ■ 
\2 


^[<l>^-4»si) 


(3.173) 


and  that 


tan^^[®  ^5  +  ®  54  j)  =  cot^-  AX, 


cos 


V2 


7C 

2 


j. 

--  +  <1>4 


cosf- 


^2 


^  j. 

+ - ^A 

2 


(3.174) 


1 


C0s(^^[<|>4  -<1>5]] 


=  C0tl-AX,  — 7- -  N 


COl 


(3.175) 


cotl  —  AX, 


COl 


(3.176) 


Equation  3.173  becomes 
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>15  “  ®  a4  ]  -  tan 


2 


co^  ^[^A  + 


(3.177) 


and  Equation  3.176  is  equivalent  to 


+^54]- tan 


2 


sinQ[<|)^  + 


(3.178) 


Therefore  we  can  find  <»  ,„  and  0)^^  fhrougji  the  following  two  formulas 


(0^5  =  tan 


cotl  —  AX 


.  cos 

^  •  1 

sin 

[\[*A  +'t'5]) 

+  tan  t-i 


1 


COSI  -[^A  + 


(3.179) 


0^4  =tan  ^<1 


^COS 

_ —  k 

sin| 

-1 


.... 

J 

COS 

Q[<}>>1  +4>s]) 

(3.180) 


Putting  these  angles  in  terms  of  heading  we  find  that  the  heading  from  point  A  to  point  B  is 


Heading  A->  B  =  (o  ab 

and  the  heading  from  B  to  A  is  the  mirror  image  of  A  to  B 

Heading  B A=^2n-(o qa 


(3.181) 


(3.182) 


Figure  3.9  illustrates  the  case  when  the  longitude  of  A  is  larger  than  the  longitude  of  B 


Ak  =  —  Xq 

AX>0;  AX  <TC 

1  ,  ,  .71 


Figure  3.9  HEADING  BETWEEN  POINTS  A  AND  B 

Here  the  same  formulas  for  ©ab  and  ©ba  apply.  The  heading  from  point  A  to  B  is 

Heading  A B  =  2%  -  a  (3.183) 

and  the  heading  from  B  to  A  is  the  mirror  image  of  A  to  B 

Heading  B  A  =  (ii  ba  (3.184) 

Now  that  we  understand  the  how  to  obtain  the  heading  between  two  points  the 
next  step  is  to  apply  this  theory  to  our  two  coordinate  systems.  Let  ©,•  denote  the  heading 

from  the  center  of  each  ellipse  <|>/,X./  to  the  north  pole^„p  in  the  ^Tji^-coordinate  system. 

To  find  ©  / ,  let  ,Xi  -  point  A  and  ^„p  =  point  B  and  execute  the  pseudo  code  in  the 

beginning  of  this  section.  If  0;  is  the  heading  of  the  ellipse  at  the  point  in  the  xyz- 

coordinate  system,  then  the  heading  of  the  ellipse  at  the  point  <|)j ,  Xj  in  the  ^T|^-coordinate 
system  is  the  sum  of  the  original  heading  and  the  offset  caused  by  rotation 

e;=®/+0/  (3.185) 
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IV.  COMBINATION  ALGORITHM 


A.  PROCESS  OUTLINE 

The  Combination  Algorithm  will  apply  the  theory  and  concepts  from  Chapters  II 
and  III  to  the  real  world  environment.  The  input  parameters  for  the  combination 
algorithm  will  consist  of  a  set  of  error  ellipses  with  center  <l>i ,  A/  (/=!,  ,  N),  semi-major 

axis  A, ,  semi-minor  axis  B, ,  heading  Oj  and  probability  pj  (=>  Xi  )  •  algorithm  outlined 
below  will  process  the  N  input  ellipses  and  produce  a  solution  ellipse  that  represents  an  area 
that  has  the  specified  probability  that  the  target  lie  within  it’s  bounds.  FoUowing  is  the 
Combination  algorithm  presented  in  a  step  by  step  manner. 

1.  Input  a  set  of  error  ellipses  with  center  ,  A,  (/=1,  ...  ,  N),  semi-major  axis  A  , 

semi-mmor  axis  B, ,  heading  Of  and  probability  (=>  ) ,  where  ,  A/  refer  to  the  latitude 

and  longitude  respectively.  The  semi-axes  A/  and  B/  are  in  an  appropriate  unit  of  length,  6i 
in  degrees  and  the  desired  probability  level  p  refers  to  the  confidence  that  the  target  lie 
within  the  generated  solution  ellipse. 

2.  Project  semi-axes  A  and  B/  onto  a  plane  to  obtain  a;  and  bj,  where 


a,  =  r  tan 


'^1 


(4.1) 


b,  =  r,  tan 


.r.J 


(4.2) 


Rationalization  and  development  of  Equation  4.1  and  Equation  4.2  is  presented  in 
Chapter  III.A.1,  Scaled  Covariance  Matrix. 

3.  Find  the  eigenvalues  corresponding  to  the  semi-major  axes  a, ,  semi-minor  axes  b, 
and  the  probability  level  p  of  each  ellipse,  where 
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(4.3) 


(4.4) 


Equation  4.3  and  Equation  4.4  find  their  bases  in  Chapter  II.A,  Relationship  of  the 
Error  Ellipse  to  a  Bivariate  Normal  Distribution.  The  values  of  p  for  the  ^  distributed 
function  are  given  in  Table  2.1. 


4.  Find  the  first  approximation  for  the  combined  location  of  the  center  of  the 


ellipses 


(4.5) 


T  =  =  ^Scos(<^)sin(^) 


(4.6) 


1  ^ 


=  ^Zsin(4 


(4.7) 


(|>  =  sin~^ 


(4.8) 


X,  =  tan  ^  — 


(4.9) 


The  first  part  of  this  step  finds  the  simple  average  of  the  (x,y,z)  components  of  the 
individual  error  ellipses.  Note  that  the  information  within  die  individual  observation 


vectors  x, ,  where  x,-  = 


T, 

\z,J 


should  be  retained  for  future  use.  These  x,  vectors  will  be 


needed  for  Step  5.  The  second  part  of  this  step  converts  the  averaged  (x,j'<//)  components, 
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producing  the  latitude  (f>  and  longitude  X  which  is  an  approximation  of  die  center  of  all  the 
observed  error  ellipses. 

5.  Transform  the  location  of  the  center  of  each  error  ellipse  from  xyz-coordinates 
to  ^74~coordinates  (prime  coordinates). 


o 

o 

o 

o 

sin(/l)  cos(^) 

sin(^)^ 

cos(^^)cos(a,P 

-  sin(A) 

cos(A) 

0 

r,  cos(<^)sin(;i,) 

(4.10) 

cos(A)  sin(^) 

-  sin(A)  sin(^) 

cos(^)j 

reSin(<^i)  ; 

^  cos(/l)  cos(^) 

sin(/l)  cos 

sin(^)^ 

-  sin(A) 

cos(/l) 

0 

yi 

(4.11) 

V-  cos(A)  sin(^) 

-  sin(A)  sin(^] 

O 

o 

U,J 

<j>i  =  sin*’ 

(4.12) 

=  tan  ’ 


Ik 


(4.13) 


Note  that  x,  calculated  in  Step  4  does  not  need  to  be  recalculated  to  obtain 


where 


Vi 


.  The  transformation  matrix  used  in  Equation  4.10  is  based  on  Equation 


3.158  developed  in  Chapter  III.C,  Ellipse  Rotation.  During  this  step  the  first  approximation 
coordinates  <f>,X  (Step  4)  are  the  focus  upon  which  the  coordinate  system  is  rotated  as 
illustrated  in  Figure  4.1.  The  end  product  being  the  spatial  relationship  of  each  observation 
<I>,X  to  the  first  approximation  <|^^,X^  is  maintained,  but  transformed  to  a  coordinate  system 
where  ^  =  0,A  =  0.  This  transformation  is  required  compensate  for  the  difference  in  arc- 
length  due  to  a  fixed  change  in  longitude  as  the  latitude  increases  firom  zero. 
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First  Approximation 
Coordinates 


6.  Find  the  coordinates  of  the  north  pole  in  ^//^-coordinates 


=  sin  '(cos^) 

(4.14) 

Jo  if  ^  >  0 

(4.15) 

~  |;r  if  <  0 

Location  of  the  north  pole  in  prime  coordinates  is  needed  to  solve  for  die  new 
ellipse  heading  after  the  coordinate  rotation.  Figure  4.2  shows  the  geometry  involved  in  the 
pole  rotation.  Equations  4.14  and  4.15  were  obtained  form  Equations  3.169  and  3.170. 
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Figure  4.2  NORTH  POLE  ROTATION 

7.  Find  the  heading  from  the  center  of  each  error  ellipse  ^ ^  to  the  north 
pole  ip  ^  ^//^-coordinates. 

The  process  for  this  step  is  detailed  in  Chapter  III.D,  Ellipse  Heading  Relative  to 
North  Pole.  In  execution,  the  first  move  is  to  examine  the  trivial  cases  ^  =  Kp  , 

^  =  A„p+7r,  and  the  case  where  one  of  the  points  is  located  on  the  north  pole.  In  the 
likely  event  that  none  of  the  trivial  case  match,  Napier’s  analogies  must  be  utilized  by 
employing  Equation  3.179. 

8..Find  the  heading  of  the  error  ellipses  in  |/7^-coordinates 

^  ^  (4-16) 

The  heading  of  the  individual  error  ellipses  in  the  prime  coordinate  system  is  the 
sum  of  the  heading  of  each  error  ellipse  (Step  7)  and  ellipse’s  heading  in  xyz-coordinates 

(Step  1). 

9.  Find  the  scaled  covariance  matrix  corresponding  to  each  error  ellipse 
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r  -  ’'M 

'  In  PJ 

(4.17) 

€-=  '  f"’-  -"1 

(4.18) 

11 

+ 

O 

o 

1 

(4.19) 

1 

O 

O 

1 

11 

(4.20) 

=sin^  cos^(£^(,)-f  (,)] 

1 

(4.21) 

Along  with  the  scaled  covariance  matrix  C,. ,  Step  10  will  require  die  individual 

ellipse  scaled  covariance  matrix  inverse  C"’  also  be  calculated.  The  elements  of  these 
matrices  are  derived  in  Chapter  III.A.3,  Variable  Relationships. 

10.  Find  the  2x2  scaled  covariance  matrix  for  the  sum  of  the  error  ellipses 


c-'=2c-  = 

i=l 


<T  V 
kO  tj 


(4.22) 


1  1 

f  r 

-u^ 

1 «  )  - 

(JT-O^ 

\-o 

cr  ) 

vr 

P) 

(4.23) 


Equation  4.22  is  a  result  of  the  methodology  outlined  in  Chapter  II.C.  The  2x2 

scaled  covariance  matrix  for  the  sum  of  the  error  elUpses  C  produced  in  this  step 
represents  the  solution  ellipse  that  has  the  specified  probability  that  the  target  lie  within  it’s 
bounds.  The  location  and  heading  of  this  solution  ellipse  has  yet  to  be  determined. 

11.  Convert  the  scaled  covariance  matrix  for  the  sum  of  the  ellipses  to  the  3x3 
matrix 
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(4.27) 


(4.28) 


There  are  foiar  tasks  to  accomplish  prior  to  producing  the  3x3  scaled  covariance 
matrix  for  the  sum  of  the  ellipses  •  First,  the  inverse  of  the  scaling  matrix  defined  in 

Equation  3.7  must  be  generated.  Second,  using  the  2x2  scaled  covariance  matrix  C  from 
Step  10,  calculate  the  unsealed  covariance  matrix  in  terms  of  latitude  and  longitude.  Next, 
populate  the  3x3  unsealed  covariance  matrix.  This  matrix  ~  is  in  terms  of  latitude, 

longitude  and  hei^t.  Lastly,  the  3x3  scaling  matrix  T  needs  to  be  produced.  T  contains 
the  partial  derivatives  of  ^,1]  and  ^  with  respect  to  X  and  h  (Equations  3.19  thru  3.27). 
With  those  tasks  complete  can  be  calculated.  This  conversion  process  is  described 

in  detail  in  Chapter  III.A.2.a  2x2  to  3x3  Covariance  Matrix  Conversion. 
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12.Find  the  wei^ted  average  location  for  the  center  of  the  error  ellipse  in 
coordinates 


M 

=  «Z«;' 

Vi 

1=1 

<f>  =  sin 


vr.J 


(4.29) 


(4.30) 


A  =  tan 


-1 


(i 

U. 


(4.31) 


The  vector  ^  represents  the  location  of  the  wei^ted  average  of  the  N 
observations.  The  theory  behind  Equation  4.29  is  presented  in  Chapter  II.C,  Methodology. 
The  new  latitude  ^  and  longitude  J  are  in  essence  refinements  of  the  rotated  first 
approximation.  The  level  of  refinement  will  be  monitored  once  the  solution  is  transformed 
back  to  the  xyz-coordinate  system.  If  convergence  upon  the  true  solution  is  not  obtained, 
an  iterative  process  will  begin. 

13.  Find  the  eigenvalues  of  the  combined  covariance  matrix. 


s,  =\[cc+p)+^a-pY  +4r"f' 

(4.32) 

(4.33) 

where  a,fi  and  y  are  obtained  from  Equation  4.23. 

The  eigenvalues  produced  in  this  step  are  the  key  to  extracting  information  about 
the  semi-major  axis,  semi-minor  axis  and  heading  of  the  covariance  matrix.  Equation  4.32 
and  Equation  4.33  supply  the  real  solutions  to  the  characteristic  polynomial  of  the 
covariance  matrix  discussed  in  Chapter  III.A.3.a,  Eigenvalues. 
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14.  Find  the  heading  of  the  error  ellipse  in  the  ^T/^^'CCordinate  system 


9 


=  tan'‘ 


S^-P) 


(4.34) 


The  error  ellipse  heading  is  derived  from  the  ^  and  7/  components  of  the 
eigenvectors  found  in  Section  III.A.3.C,  Eigenvectors.  The  development  of  Ecjuation  4.34 
can  be  found  in  Chapter  III.A.4,  Variable  Relationships. 

15.  Find  the  heading  m  from  (!> ,  X  to  X^. 


The  same  process  found  in  Step  7  is  utilized  here. 

16.  Find  the  heading  of  the  error  ellipse  in  the  xyz-coordinate  system 

9=9  -  S  (4.35) 


The  heading  of  the  error  ellipse  in  the  xyz-coordinate  system  is  found  subtracting 
the  headbg  offset  due  to  rotation  from  the  ellipses  heading  in  the  prime  coordinate  system. 
Note  that  Equation  4.35  agrees  with  Equation  4.16  from  Step  8. 

17.  Find  the  coordinates  of  the  wei^ted  average  location  for  the  center  of  the 
error  ellipses  in  xyz-coordinates 


fi) 

^cos(A)cos(^) 

-  sin(A) 

1 

O 

o 

y 

= 

sin(A)cos(^) 

cos(A) 

-sin(yl)sin(^^) 

U; 

1  sin(<^) 

0 

o 

o 

Id 

=  sin 


-t 


T 

UeJ 


(3.36) 


(4.37) 


A  =  tan  ‘ 


(4.38) 
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The  transformation  from  prime  coordinates  to  the  xyz-coordinate  system  is 
developed  in  Chapter  III.C,  ElUpse  Rotation.  This  step  is  the  inverse  of  Step  5. 

18.  Check  for  convergence 


ff  ^  I  <  and 

else  set 


X-X 


<K]^  then 


goto  Step  19 

X  =  X  ;  y=y  ;  z=z 


(f>  '=■  <f>  X  —  X 

and  goto  Step  5 


In  some  cases  the  difference  between  the  first  approximation  and  the 
weighted  average  location  my  be  significant.  By  setting  values  for  and  that  are  small 


an  iterative  process  is  established  to  fine  tune  the  solution. 

19.  Calculate  the  semi-major  and  semi-minor  axes  of  the  combined  error  ellipse  for 
the  derived  probability  level 


a  = 

b=>/?r 


(4.39) 

(4.40) 


The  semi-major  and  semi-minor  axes  of  the  combined  error  ellipse  are  found  by 
utilizing  the  eigenvalues  calculated  in  Step  14  and  a  selected  chi-square  number  for  Table 
II.1.  Chapter  II.A,  Relationship  of  the  Error  ElUpse  to  A  Bivariate  Normal  Distribution 
gives  the  origins  of  Equation  4.39  and  Equation  4.40. 

20.  Project  the  semi-major  and  semi-minor  axes  from  the  plane  to  a  spherical  earth 


A  = 


r,tan' 


(4.41) 
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The  last  step  in  the  Combination  Algorithm  places  solution  ellipse  back  onto  the 
spherical  earth.  This  reverses  Step  2  and  supplies  the  last  parameters  of  the  solution  ellipse 

which  consists  of  latitude  ^ ,  longitude  X ,  semi-major  axes  A,  semi-minor  axes  B,  and 

heading  9. 

B.  PROCESS  FLOW  CHART 

Figure  4.3  provides  a  condensed  view  of  the  algorithm.  For  details  of  specific  steps 
refer  to  Section  A  of  this  chapter. 


Figute  4.3  PROCESS  FLOW  CHART 


62 


V.  SUMMARY  &  CONCLUSIONS 


A.  AREAS  FOR  FUTURE  STUDY 

One  of  the  key  assumptions  made  in  the  development  of  this  algorithm  was  that  the 
error  ellipse  data  is  derived  from  a  bivariate  normal  distribution.  Several  transformations 
were  made  throu^out  the  derivation  to  transform  spherical  data  to  a  form  that  could  be 
analyzed  using  a  bivariate  normal  distribution.  When  the  error  ellipse  is  small  (a  «  r_e)  the 
probability  distribution  is  hi^ly  localized  and  these  transformations  lead  to  a  good 
approximation  of  a  spherical  probability  distribution.  The  algorithm  developed  here  will 
therefore  be  applicable  to  most  error  ellipses  obtained  from  operational  systems.  Strictly 
speaking,  however,  we  should  have  used  a  probability  distribution  more  suited  to  spherical 
data,  the  Kent  distribution  (Fisher,  Lewis,  Embletonl987).  While  this  distribution  is  more 
difficult  to  work  with,  it  would  be  valid  for  error  ellipses  of  all  sizes.  This  issue  should  be 
examined  more  fully  in  the  future. 

There  are  three  remaining  areas  relevant  to  this  thesis  that  are  available  for  future 
research.  The  first  and  foremost  is  to  design  and  implement  a  test  and  validation  program. 
The  validation  program  would  provide  a  means  to  quantitatively  verify  the  theory  put  forth 
in  this  research.  The  project  would  entail  producing  computer  code  to  implement  the 
Combination  Algorithm.  Real  world  data  sets  are  available  to  compare  against  output  that 
would  be  generated  by  the  program. 

Broadening  the  Combination  Algorithm  to  include  multiple  source  data  is  the  next 
potential  research  area.  The  ability  to  correlate  observations  from  multiple  platforms  would 
greatly  enhance  our  national  warfi^ting  capability.  Decision  makers  would  no  longer  need 
to  mentally  extrapolate  a  solution  from  numerous  observations  of  the  same  target.  The 
generation  of  this  composite  geolocation  would  reduce  total  decision  cycle  time,  ultimately 
leading  to  increased  force  effectiveness.  With  an  effective  combination  algorithm  m  place 
the  need  to  transmit  multiple  geolocations  of  the  same  target  would  be  eliminated,  reducing 
fleet  bandwidth  requirements.  By  combining  observations  of  the  various  platforms  both 
bandwidth  and  data  storage  requirements  are  reduced.  The  reduction  would  be  a  result  of 
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two  factors:  (1)  collection  platforms  will  combine  same  source  observations  reducing  the 
amount  of  data  to  be  transmitted,  (2)  once  the  collection  platform  solutions  are  received 
and  combined  into  a  composite  solution  the  data  received  from  the  collection  platforms  can 
be  discarded.  The  problem  of  including  multiple  source  data  could  be  approached  by 
expanding  the  observation  vector  to  include  system  dependent  biases  such  as  atmospheric 

delays  and  refraction. 

The  third  area  to  consider  for  future  research  is  the  correlation  of  observations 
obtained  of  a  moving  target.  Expanding  the  domain  of  possible  targets  from  fixed  site  to 
include  mobile  platforms  would  significantly  reduce  the  information  presented  to  a  decision 
maker.  By  expanding  the  number  of  platforms  available  for  data  reduction  through 
combination,  bandwidth  and  data  storage  resources  are  conserved.  The  inclusion  of 
moving  targets  might  be  accompUshed  by  sampling  and  processing  a  discrete  number  of  hits 
over  a  large  number  of  observations.  Those  interested  in  pursuing  any  one  of  these  area 
can  contact  NRaD  Code  841  for  information  on  current  research. 

Areas  of  interest  observed  during  initial  attempts  to  implement  the  Combination 
Algorithm  revolved  around  data  structure  design.  Because  the  vast  majority  of  operations 
to  be  performed  in  the  algorithm  involve  matrix  arithmetic,  a  language  designed  to  for  that 
purpose  such  as  MATLAB  is  recommended.  If  C++  is  to  be  used,  it  is  recommended  that 
an  estabUshed  toolbox  of  data  objects  and  algorithms  such  as  M++  by  Dayd  Software  be 
utilized.  Because  of  the  large  number  of  multi-dimensioned  arrays  required  to  code  this 
algorithm  a  large  memory  model  will  be  required.  In  the  event  that  a  smaller  memory 
model  must  be  used,  the  number  of  observations  the  program  wUl  be  able  to  process  will 
be  severely  restricted.  To  further  reduce  the  strain  on  memory  resources  dynamic  allocation 
of  objects  should  be  incorporated. 
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B.  CONCLUSION 

This  thesis  has  set  the  theoretical  basis  and  developed  an  algorithm  for  optimally 
combining  multiple  small  spherical  error  ellipses  into  a  single  error  ellipse.  This  algorithm  is 
appUcable  to  most  of  the  error  ellipses  produced  by  operational  systems.  When  this 
algorithm  is  validated  it  could  have  a  major  impact  on  how  geolocation  data  is  presented  to 
the  warfigjiter.  In  the  future,  as  the  electronic  density  of  the  battle  field  increases,  the  need 
for  fusion  of  geolocation  data  will  become  increasin^y  apparent.  This  work  can  be 
considered  the  first  step  toward  providing  die  warfi^ter  with  a  fused  multi-sensor 
geolocation  solution. 
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